% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Replication "Deconstructing the Yield Curve"
% Crump and Gospodinov (2024)
% Date: 26-JUL-2024
% Code to produce Table 9
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all; 
clc; 
close all;
rng(12345);
addpath('functions');
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Options
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
o.exactDataBR           = false; 
o.dataSet               = 'gsw';
o.freq                  = 'q';
o.startDateSample       = 197104; 
o.endDateSample         = 201801;
o.rstarBasedOnCPI       = true;
p.N                     = 60;
p.yldMats               = [1,2,4:4:60];
p.H                     = 1;
p.CG_B_BC               = 1000;
p.BH_B_BC               = 1000;
p.B                     = 5000;
p.L                     = 3;
p.M                     = 'ROT';
p.K                     = numel(p.yldMats);
% Keep this fixed to ensure we have yoy core PCE starting in 1961Q4 for EWMA
p.startDate             = 196004; 
idxSaveData             = true;
nwLagh4                 = 6;

if o.exactDataBR && (o.startDateSample ~= 197104 ||  o.endDateSample > 201801)
    error('Option to use exact BR data only for sample: 197104--201801');
end

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Load data
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp(' '); disp('Loading Data...'); disp(' ');
% Load data sets
dataAll = loadBondData(p.N);
dataTag = ['dataAll.', o.dataSet, '.', o.freq];
for i = {'prc', 'yld', 'ret', 'xret', 'fwd', 'dret'}, eval([i{1}, 'Raw = ', dataTag, '.', i{1}, ';']);  end
data.mats = (1:p.N);
datesBonds = eval([dataTag, '.dates']);
data.period = eval([dataTag, '.period']);

datesAll = 100*kron((1901:2026)',ones(4,1))+kron(ones(126,1),(1:4)'); %!%
tmpIdxStart = find(datesAll==p.startDate);
tmpIdxEnd = find(datesAll==o.endDateSample);
datesAll = datesAll(tmpIdxStart:tmpIdxEnd);
TT = numel(datesAll);

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Load CPI data
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
rawCPI = readtable('../input/CPILFESL.xls', 'DataRange', 'A11');
% Lag CPI by one month to match Cieslak and Pavola approach
rawCPI.CPILFESL = [NaN; rawCPI.CPILFESL(1:end-1)];
rawCPI.coreCPIyoy = [NaN(12,1); 100*(log(rawCPI.CPILFESL(13:end))-log(rawCPI.CPILFESL(1:end-12)))];
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Load PTR data
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
rawPTR = readtable('../input/115622-V1/data/pistar_PTR.csv');
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Load PCE data
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
rawPCE = readtable('../input/115622-V1/data/PCEPILFE.csv');
rawPCE.corePCEyoy = [NaN(12,1); 100*(log(rawPCE.PCEPILFE(13:end))-log(rawPCE.PCEPILFE(1:end-12)))];
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Load T-Bill data
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if o.exactDataBR
    rawTB3mavg = readtable('../input/115622-V1/data/TB3MS.csv');
end
rawTB3 = readtable('../input/DTB3.xls', 'Range', 'A11');
rawTB6 = readtable('../input/DTB6.xls', 'Range', 'A11');

% PTR
rawDatesPTR = 100*rawPTR.Year + rawPTR.Quarter;
data.PTR = matchDates(datesAll,rawDatesPTR,rawPTR.pistar_PTR);
% PCE
tmpIdx = find(month(rawPCE.DATE)==3, 1, 'first');
tmpData = rawPCE.corePCEyoy(tmpIdx:3:end);
rawDatesPCE = 100*year(rawPCE.DATE(tmpIdx:3:end)) + quarter(rawPCE.DATE(tmpIdx:3:end));
data.corePCEyoy = matchDates(datesAll,rawDatesPCE,tmpData);
% CPI
tmpIdx = find(month(rawCPI.observation_date)==3, 1, 'first');
tmpData = rawCPI.coreCPIyoy(tmpIdx:3:end);
rawDatesCPI = 100*year(rawCPI.observation_date(tmpIdx:3:end)) + quarter(rawCPI.observation_date(tmpIdx:3:end));
data.coreCPIyoy = matchDates(datesAll,rawDatesCPI,tmpData);
% TB3
tmpIdx = find(month(rawTB3.observation_date)==3, 1, 'first');
tmpData = rawTB3.DTB3(tmpIdx:3:end);
rawDatesTB3 = 100*year(rawTB3.observation_date(tmpIdx:3:end)) + quarter(rawTB3.observation_date(tmpIdx:3:end));
data.TB3 = matchDates(datesAll,rawDatesTB3,tmpData);
if o.exactDataBR
    tmpIdx = find(month(rawTB3mavg.DATE)==3, 1, 'first');
    tmpData = rawTB3mavg.TB3MS(tmpIdx:3:end);
    rawDatesTB3mavg = 100*year(rawTB3mavg.DATE(tmpIdx:3:end)) + quarter(rawTB3mavg.DATE(tmpIdx:3:end));
    data.TB3mavg = matchDates(datesAll,rawDatesTB3mavg,tmpData);
end
% TB6
tmpIdx = find(month(rawTB6.observation_date)==3, 1, 'first');
tmpData = rawTB6.DTB6(tmpIdx:3:end);
rawDatesTB6 = 100*year(rawTB6.observation_date(tmpIdx:3:end)) + quarter(rawTB6.observation_date(tmpIdx:3:end));
data.TB6 = matchDates(datesAll,rawDatesTB6,tmpData);
% Bond Data (replace 3-month and 6-month with Fred data)
data.yld = matchDates(datesAll,datesBonds,yldRaw(:,p.yldMats));
data.yld(:,1) = data.TB3;
data.yld(:,2) = data.TB6;

if o.exactDataBR
    warning('OFF', 'MATLAB:table:ModifiedAndSavedVarnames');    
    rawY = readtable('../input/115622-V1/data/yields.csv'); 
    yearY = floor(rawY.Date/10000);
    monthY = floor((rawY.Date-10000*floor(rawY.Date/10000))/100);
    quarterY = monthY/3;
    datesY = 100*yearY+quarterY;
    tmpY = table2array(rawY(:,2:end));
    if datesY(end)>datesAll(end)
        tmpIdx = find(datesY==datesAll(end));
        datesY = datesY(1:tmpIdx);
        tmpY = tmpY(1:tmpIdx,:);
    end
    [~,tmpIdx] = intersect(datesAll,datesY);        
    tmpYupdate = data.yld(tmpIdx,:);
    tmpYupdate(~isnan(tmpY)) = tmpY(~isnan(tmpY));
    data.yld(tmpIdx,:) = tmpYupdate;
end

[xrn,xrn4] = makeReturnsBR(data.yld);

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Make PCs
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
idxSampleStart = find(datesAll==o.startDateSample);
idxSampleEnd = find(datesAll==o.endDateSample);
T = idxSampleEnd - idxSampleStart + 1;
% Make PCs
[tmp_evecs,tmp_eigs] = eig(cov(data.yld(idxSampleStart:idxSampleEnd,:)));
[~,idx] = sort(diag(tmp_eigs), 'descend');
tmp_evecs = tmp_evecs(:,idx);
PCsLoad = tmp_evecs(:,1:3)*diag(sign(tmp_evecs(end,1:3)));
PCs = data.yld(idxSampleStart:idxSampleEnd,:)*PCsLoad;
PC1 = PCs(:,1);
PC2 = PCs(:,2);
PC3 = PCs(:,3);
% Make scaled PCs
sc = [sum(tmp_evecs(:,1)), tmp_evecs(1,2)-tmp_evecs(end,2), tmp_evecs(end,3)-2*tmp_evecs(4,3)+tmp_evecs(1,3)];
PCssc = PCs*diag(1./sc);
PC1sc = PCssc(:,1);
PC2sc = PCssc(:,2);
PC3sc = PCssc(:,3);
ehat = data.yld(idxSampleStart:idxSampleEnd,:) - PCs*PCsLoad';
s2hat = trace(ehat'*ehat)/numel(ehat);

xrnSample = xrn(idxSampleStart:idxSampleEnd,4:end);
xrn4Sample = xrn4(idxSampleStart:idxSampleEnd,4:end);
PTRSample = data.PTR(idxSampleStart:idxSampleEnd);

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Make all other variables
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if ~o.rstarBasedOnCPI
    if  o.exactDataBR
        RRQ = data.TB3mavg - data.corePCEyoy;
    else
        RRQ = data.TB3 - data.corePCEyoy;
    end
else
    RRQ = data.TB3 - data.coreCPIyoy;
end
ewmaRRQ = [NaN(4,1); ewma(RRQ(5:end),0.98)];
ewmaRRQsample = ewmaRRQ(idxSampleStart:idxSampleEnd);

CPoQ = ewma_trunc(data.coreCPIyoy,0.9615,40);
CPoQsample = CPoQ(idxSampleStart:idxSampleEnd);

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% All regression results
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
betaHat1 = NaN(5,15);
betaHat2 = NaN(6,15);
betaHat1_CPo = NaN(5,15);
betaHat2_CPo = NaN(6,15);
betaHat3 = NaN(5,15);
betaHat41 = NaN(5,15);
betaHat42 = NaN(6,15);
betaHat41_CPo = NaN(5,15);
betaHat42_CPo = NaN(6,15);
betaHat43 = NaN(5,15);
%
tStat1 = NaN(5,15);
tStat2 = NaN(6,15);
tStat1_CPo = NaN(5,15);
tStat2_CPo = NaN(6,15);
tStat3 = NaN(5,15);
tStat41 = NaN(5,15);
tStat42 = NaN(6,15);
tStat41_CPo = NaN(5,15);
tStat42_CPo = NaN(6,15);
tStat43 = NaN(5,15);
%
R2_1 = NaN(15,1);
R2_1_CPo = NaN(15,1);
R2_2 = NaN(15,1);
R2_2_CPo = NaN(15,1);
R2_3 = NaN(15,1);
R2_40 = NaN(15,1);
R2_41 = NaN(15,1);
R2_41_CPo = NaN(15,1);
R2_42 = NaN(15,1);
R2_42_CPo = NaN(15,1);
R2_43 = NaN(15,1);
%
betaHat1_sub = NaN(5,15);
betaHat2_sub = NaN(6,15);
betaHat1_sub_CPo = NaN(5,15);
betaHat2_sub_CPo = NaN(6,15);
betaHat3_sub = NaN(5,15);
betaHat41_sub = NaN(5,15);
betaHat42_sub = NaN(6,15);
betaHat41_sub_CPo = NaN(5,15);
betaHat42_sub_CPo = NaN(6,15);
betaHat43_sub = NaN(5,15);
%
tStat1_sub = NaN(5,15);
tStat2_sub = NaN(6,15);
tStat1_sub_CPo = NaN(5,15);
tStat2_sub_CPo = NaN(6,15);
tStat3_sub = NaN(5,15);
tStat41_sub = NaN(5,15);
tStat42_sub = NaN(6,15);
tStat41_sub_CPo = NaN(5,15);
tStat42_sub_CPo = NaN(6,15);
tStat43_sub = NaN(5,15);
%
R2_1_sub = NaN(15,1);
R2_1_sub_CPo = NaN(15,1);
R2_2_sub = NaN(15,1);
R2_2_sub_CPo = NaN(15,1);
R2_3_sub = NaN(15,1);
R2_40_sub = NaN(15,1);
R2_41_sub = NaN(15,1);
R2_41_sub_CPo = NaN(15,1);
R2_42_sub = NaN(15,1);
R2_42_sub_CPo = NaN(15,1);
R2_43_sub = NaN(15,1);

offset = find(datesAll==198501) - idxSampleStart;

for i = 1:14
    
    reg1 = nwest(xrnSample(2:end,i),[ones(T-1,1), PC1sc(1:end-1,:), PC2sc(1:end-1,:), PC3sc(1:end-1,:) PTRSample(1:end-1)],0);
    reg1_CPo = nwest(xrnSample(2:end,i),[ones(T-1,1), PC1sc(1:end-1,:), PC2sc(1:end-1,:), PC3sc(1:end-1,:) CPoQsample(1:end-1)],0);
    reg41 = nwest(xrn4Sample(5:end,i),[ones(T-4,1), PC1sc(1:end-4,:), PC2sc(1:end-4,:), PC3sc(1:end-4,:) PTRSample(1:end-4)],nwLagh4); 
    reg41_CPo = nwest(xrn4Sample(5:end,i),[ones(T-4,1), PC1sc(1:end-4,:), PC2sc(1:end-4,:), PC3sc(1:end-4,:) CPoQsample(1:end-4)],nwLagh4);     
    reg2 = nwest(xrnSample(2:end,i),[ones(T-1,1), PC1sc(1:end-1,:), PC2sc(1:end-1,:), PC3sc(1:end-1,:) PTRSample(1:end-1) ewmaRRQsample(1:end-1)],0);
    reg2_CPo = nwest(xrnSample(2:end,i),[ones(T-1,1), PC1sc(1:end-1,:), PC2sc(1:end-1,:), PC3sc(1:end-1,:) CPoQsample(1:end-1) ewmaRRQsample(1:end-1)],0);
    reg42 = nwest(xrn4Sample(5:end,i),[ones(T-4,1), PC1sc(1:end-4,:), PC2sc(1:end-4,:), PC3sc(1:end-4,:) PTRSample(1:end-4) ewmaRRQsample(1:end-4)],nwLagh4); 
    reg42_CPo = nwest(xrn4Sample(5:end,i),[ones(T-4,1), PC1sc(1:end-4,:), PC2sc(1:end-4,:), PC3sc(1:end-4,:) CPoQsample(1:end-4) ewmaRRQsample(1:end-4)],nwLagh4);     
    reg3 = nwest(xrnSample(2:end,i),[ones(T-1,1), PC1sc(1:end-1,:), PC2sc(1:end-1,:), PC3sc(1:end-1,:) ewmaRRQsample(1:end-1)],0);
    reg43 = nwest(xrn4Sample(5:end,i),[ones(T-4,1), PC1sc(1:end-4,:), PC2sc(1:end-4,:), PC3sc(1:end-4,:) ewmaRRQsample(1:end-4)],nwLagh4); 
    %
    betaHat1(:,i) = reg1.beta;
    betaHat1_CPo(:,i) = reg1_CPo.beta;
    betaHat2(:,i) = reg2.beta;
    betaHat2_CPo(:,i) = reg2_CPo.beta;
    betaHat3(:,i) = reg3.beta;
    betaHat41(:,i) = reg41.beta;
    betaHat41_CPo(:,i) = reg41_CPo.beta;
    betaHat42(:,i) = reg42.beta;
    betaHat42_CPo(:,i) = reg42_CPo.beta;
    betaHat43(:,i) = reg43.beta;
    %
    tStat1(:,i) = reg1.tstat;
    tStat1_CPo(:,i) = reg1_CPo.tstat;
    tStat2(:,i) = reg2.tstat;
    tStat2_CPo(:,i) = reg2_CPo.tstat;
    tStat3(:,i) = reg3.tstat;
    tStat41(:,i) = reg41.tstat;
    tStat41_CPo(:,i) = reg41_CPo.tstat;
    tStat42(:,i) = reg42.tstat;
    tStat42_CPo(:,i) = reg42_CPo.tstat;
    tStat43(:,i) = reg43.tstat;   
    %
    R2_1(i) = reg1.rbar;
    R2_1_CPo(i) = reg1_CPo.rbar;
    R2_2(i) = reg2.rbar;
    R2_2_CPo(i) = reg2_CPo.rbar;
    R2_3(i) = reg3.rbar;
    R2_41(i) = reg41.rbar;
    R2_41_CPo(i) = reg41_CPo.rbar;
    R2_42(i) = reg42.rbar;
    R2_42_CPo(i) = reg42_CPo.rbar;
    R2_43(i) = reg43.rbar;   
    %
    clearvars reg*
    %
    reg1_sub = nwest(xrnSample(2+offset:end,i),[ones(T-1-offset,1), PC1sc(1+offset:end-1,:), PC2sc(1+offset:end-1,:), PC3sc(1+offset:end-1,:) PTRSample(1+offset:end-1)],0);
    reg1_sub_CPo = nwest(xrnSample(2+offset:end,i),[ones(T-1-offset,1), PC1sc(1+offset:end-1,:), PC2sc(1+offset:end-1,:), PC3sc(1+offset:end-1,:) CPoQsample(1+offset:end-1)],0);
    reg41_sub = nwest(xrn4Sample(5+offset:end,i),[ones(T-4-offset,1), PC1sc(1+offset:end-4,:), PC2sc(1+offset:end-4,:), PC3sc(1+offset:end-4,:) PTRSample(1+offset:end-4)],nwLagh4); 
    reg41_sub_CPo = nwest(xrn4Sample(5+offset:end,i),[ones(T-4-offset,1), PC1sc(1+offset:end-4,:), PC2sc(1+offset:end-4,:), PC3sc(1+offset:end-4,:) CPoQsample(1+offset:end-4)],nwLagh4);     
    reg2_sub = nwest(xrnSample(2+offset:end,i),[ones(T-1-offset,1), PC1sc(1+offset:end-1,:), PC2sc(1+offset:end-1,:), PC3sc(1+offset:end-1,:) PTRSample(1+offset:end-1) ewmaRRQsample(1+offset:end-1)],0);
    reg2_sub_CPo = nwest(xrnSample(2+offset:end,i),[ones(T-1-offset,1), PC1sc(1+offset:end-1,:), PC2sc(1+offset:end-1,:), PC3sc(1+offset:end-1,:) CPoQsample(1+offset:end-1) ewmaRRQsample(1+offset:end-1)],0);
    reg42_sub = nwest(xrn4Sample(5+offset:end,i),[ones(T-4-offset,1), PC1sc(1+offset:end-4,:), PC2sc(1+offset:end-4,:), PC3sc(1+offset:end-4,:) PTRSample(1+offset:end-4) ewmaRRQsample(1+offset:end-4)],nwLagh4); 
    reg42_sub_CPo = nwest(xrn4Sample(5+offset:end,i),[ones(T-4-offset,1), PC1sc(1+offset:end-4,:), PC2sc(1+offset:end-4,:), PC3sc(1+offset:end-4,:) CPoQsample(1+offset:end-4) ewmaRRQsample(1+offset:end-4)],nwLagh4);     
    reg3_sub = nwest(xrnSample(2+offset:end,i),[ones(T-1-offset,1), PC1sc(1+offset:end-1,:), PC2sc(1+offset:end-1,:), PC3sc(1+offset:end-1,:) ewmaRRQsample(1+offset:end-1)],0);
    reg43_sub = nwest(xrn4Sample(5+offset:end,i),[ones(T-4-offset,1), PC1sc(1+offset:end-4,:), PC2sc(1+offset:end-4,:), PC3sc(1+offset:end-4,:) ewmaRRQsample(1+offset:end-4)],nwLagh4); 
    %
    betaHat1_sub(:,i) = reg1_sub.beta;
    betaHat1_sub_CPo(:,i) = reg1_sub_CPo.beta;
    betaHat2_sub(:,i) = reg2_sub.beta;
    betaHat2_sub_CPo(:,i) = reg2_sub_CPo.beta;
    betaHat3_sub(:,i) = reg3_sub.beta;
    betaHat41_sub(:,i) = reg41_sub.beta;
    betaHat41_sub_CPo(:,i) = reg41_sub_CPo.beta;
    betaHat42_sub(:,i) = reg42_sub.beta;
    betaHat42_sub_CPo(:,i) = reg42_sub_CPo.beta;
    betaHat43_sub(:,i) = reg43_sub.beta;
    %
    tStat1_sub(:,i) = reg1_sub.tstat;
    tStat1_sub_CPo(:,i) = reg1_sub_CPo.tstat;
    tStat2_sub(:,i) = reg2_sub.tstat;
    tStat2_sub_CPo(:,i) = reg2_sub_CPo.tstat;
    tStat3_sub(:,i) = reg3_sub.tstat;
    tStat41_sub(:,i) = reg41_sub.tstat;
    tStat41_sub_CPo(:,i) = reg41_sub_CPo.tstat;
    tStat42_sub(:,i) = reg42_sub.tstat;
    tStat42_sub_CPo(:,i) = reg42_sub_CPo.tstat;
    tStat43_sub(:,i) = reg43_sub.tstat;   
    %
    R2_1_sub(i) = reg1_sub.rbar;
    R2_1_sub_CPo(i) = reg1_sub_CPo.rbar;
    R2_2_sub(i) = reg2_sub.rbar;
    R2_2_sub_CPo(i) = reg2_sub_CPo.rbar;
    R2_3_sub(i) = reg3_sub.rbar;
    R2_41_sub(i) = reg41_sub.rbar;
    R2_41_sub_CPo(i) = reg41_sub_CPo.rbar;
    R2_42_sub(i) = reg42_sub.rbar;
    R2_42_sub_CPo(i) = reg42_sub_CPo.rbar;
    R2_43_sub(i) = reg43_sub.rbar;       
    %
    clearvars reg*
    
end

reg1 = nwest(mean(xrnSample(2:end,:),2),[ones(T-1,1), PC1sc(1:end-1,:), PC2sc(1:end-1,:), PC3sc(1:end-1,:) PTRSample(1:end-1)],0);
reg1_CPo = nwest(mean(xrnSample(2:end,:),2),[ones(T-1,1), PC1sc(1:end-1,:), PC2sc(1:end-1,:), PC3sc(1:end-1,:) CPoQsample(1:end-1)],0);
reg41 = nwest(mean(xrn4Sample(5:end,:),2),[ones(T-4,1), PC1sc(1:end-4,:), PC2sc(1:end-4,:), PC3sc(1:end-4,:) PTRSample(1:end-4)],nwLagh4); 
reg41_CPo = nwest(mean(xrn4Sample(5:end,:),2),[ones(T-4,1), PC1sc(1:end-4,:), PC2sc(1:end-4,:), PC3sc(1:end-4,:) CPoQsample(1:end-4)],nwLagh4); 

reg2 = nwest(mean(xrnSample(2:end,:),2),[ones(T-1,1), PC1sc(1:end-1,:), PC2sc(1:end-1,:), PC3sc(1:end-1,:) PTRSample(1:end-1) ewmaRRQsample(1:end-1)],0);
reg2_CPo = nwest(mean(xrnSample(2:end,:),2),[ones(T-1,1), PC1sc(1:end-1,:), PC2sc(1:end-1,:), PC3sc(1:end-1,:) CPoQsample(1:end-1) ewmaRRQsample(1:end-1)],0);
reg42 = nwest(mean(xrn4Sample(5:end,:),2),[ones(T-4,1), PC1sc(1:end-4,:), PC2sc(1:end-4,:), PC3sc(1:end-4,:) PTRSample(1:end-4) ewmaRRQsample(1:end-4)],nwLagh4); 
reg42_CPo = nwest(mean(xrn4Sample(5:end,:),2),[ones(T-4,1), PC1sc(1:end-4,:), PC2sc(1:end-4,:), PC3sc(1:end-4,:) CPoQsample(1:end-4) ewmaRRQsample(1:end-4)],nwLagh4); 

reg3 = nwest(mean(xrnSample(2:end,:),2),[ones(T-1,1), PC1sc(1:end-1,:), PC2sc(1:end-1,:), PC3sc(1:end-1,:) ewmaRRQsample(1:end-1)],0);
reg43 = nwest(mean(xrn4Sample(5:end,:),2),[ones(T-4,1), PC1sc(1:end-4,:), PC2sc(1:end-4,:), PC3sc(1:end-4,:) ewmaRRQsample(1:end-4)],nwLagh4); 

betaHat1(:,end) = reg1.beta;
betaHat1_CPo(:,end) = reg1_CPo.beta;
betaHat2(:,end) = reg2.beta;
betaHat2_CPo(:,end) = reg2_CPo.beta;
betaHat3(:,end) = reg3.beta;
betaHat41(:,end) = reg41.beta;
betaHat41_CPo(:,end) = reg41_CPo.beta;
betaHat42(:,end) = reg42.beta;
betaHat42_CPo(:,end) = reg42_CPo.beta;
betaHat43(:,end) = reg43.beta;
%
tStat1(:,end) = reg1.tstat;
tStat1_CPo(:,end) = reg1_CPo.tstat;
tStat2(:,end) = reg2.tstat;
tStat2_CPo(:,end) = reg2_CPo.tstat;
tStat3(:,end) = reg3.tstat;
tStat41(:,end) = reg41.tstat;
tStat41_CPo(:,end) = reg41_CPo.tstat;
tStat42(:,end) = reg42.tstat;
tStat42_CPo(:,end) = reg42_CPo.tstat;
tStat43(:,end) = reg43.tstat;    
%
R2_1(end) = reg1.rbar;
R2_1_CPo(end) = reg1_CPo.rbar;
R2_2(end) = reg2.rbar;
R2_2_CPo(end) = reg2_CPo.rbar;
R2_3(end) = reg3.rbar;
R2_41(end) = reg41.rbar;
R2_41_CPo(end) = reg41_CPo.rbar;
R2_42(end) = reg42.rbar;
R2_42_CPo(end) = reg42_CPo.rbar;
R2_43(end) = reg43.rbar;
%
clearvars reg*
%
reg1_sub = nwest(mean(xrnSample(2+offset:end,:),2),[ones(T-1-offset,1), PC1sc(1+offset:end-1,:), PC2sc(1+offset:end-1,:), PC3sc(1+offset:end-1,:) PTRSample(1+offset:end-1)],0);
reg1_sub_CPo = nwest(mean(xrnSample(2+offset:end,:),2),[ones(T-1-offset,1), PC1sc(1+offset:end-1,:), PC2sc(1+offset:end-1,:), PC3sc(1+offset:end-1,:) CPoQsample(1+offset:end-1)],0);
reg41_sub = nwest(mean(xrn4Sample(5+offset:end,:),2),[ones(T-4-offset,1), PC1sc(1+offset:end-4,:), PC2sc(1+offset:end-4,:), PC3sc(1+offset:end-4,:) PTRSample(1+offset:end-4)],nwLagh4); 
reg41_sub_CPo = nwest(mean(xrn4Sample(5+offset:end,:),2),[ones(T-4-offset,1), PC1sc(1+offset:end-4,:), PC2sc(1+offset:end-4,:), PC3sc(1+offset:end-4,:) CPoQsample(1+offset:end-4)],nwLagh4); 
reg2_sub = nwest(mean(xrnSample(2+offset:end,:),2),[ones(T-1-offset,1), PC1sc(1+offset:end-1,:), PC2sc(1+offset:end-1,:), PC3sc(1+offset:end-1,:) PTRSample(1+offset:end-1) ewmaRRQsample(1+offset:end-1)],0);
reg2_sub_CPo = nwest(mean(xrnSample(2+offset:end,:),2),[ones(T-1-offset,1), PC1sc(1+offset:end-1,:), PC2sc(1+offset:end-1,:), PC3sc(1+offset:end-1,:) CPoQsample(1+offset:end-1) ewmaRRQsample(1+offset:end-1)],0);
reg42_sub = nwest(mean(xrn4Sample(5+offset:end,:),2),[ones(T-4-offset,1), PC1sc(1+offset:end-4,:), PC2sc(1+offset:end-4,:), PC3sc(1+offset:end-4,:) PTRSample(1+offset:end-4) ewmaRRQsample(1+offset:end-4)],nwLagh4); 
reg42_sub_CPo = nwest(mean(xrn4Sample(5+offset:end,:),2),[ones(T-4-offset,1), PC1sc(1+offset:end-4,:), PC2sc(1+offset:end-4,:), PC3sc(1+offset:end-4,:) CPoQsample(1+offset:end-4) ewmaRRQsample(1+offset:end-4)],nwLagh4); 
reg3_sub = nwest(mean(xrnSample(2+offset:end,:),2),[ones(T-1-offset,1), PC1sc(1+offset:end-1,:), PC2sc(1+offset:end-1,:), PC3sc(1+offset:end-1,:) ewmaRRQsample(1+offset:end-1)],0);
reg43_sub = nwest(mean(xrn4Sample(5+offset:end,:),2),[ones(T-4-offset,1), PC1sc(1+offset:end-4,:), PC2sc(1+offset:end-4,:), PC3sc(1+offset:end-4,:) ewmaRRQsample(1+offset:end-4)],nwLagh4); 
%
betaHat1_sub(:,end) = reg1_sub.beta;
betaHat1_sub_CPo(:,end) = reg1_sub_CPo.beta;
betaHat2_sub(:,end) = reg2_sub.beta;
betaHat2_sub_CPo(:,end) = reg2_sub_CPo.beta;
betaHat3_sub(:,end) = reg3_sub.beta;
betaHat41_sub(:,end) = reg41_sub.beta;
betaHat41_sub_CPo(:,end) = reg41_sub_CPo.beta;
betaHat42_sub(:,end) = reg42_sub.beta;
betaHat42_sub_CPo(:,end) = reg42_sub_CPo.beta;
betaHat43_sub(:,end) = reg43_sub.beta;
%
tStat1_sub(:,end) = reg1_sub.tstat;
tStat1_sub_CPo(:,end) = reg1_sub_CPo.tstat;
tStat2_sub(:,end) = reg2_sub.tstat;
tStat2_sub_CPo(:,end) = reg2_sub_CPo.tstat;
tStat3_sub(:,end) = reg3_sub.tstat;
tStat41_sub(:,end) = reg41_sub.tstat;
tStat41_sub_CPo(:,end) = reg41_sub_CPo.tstat;
tStat42_sub(:,end) = reg42_sub.tstat;
tStat42_sub_CPo(:,end) = reg42_sub_CPo.tstat;
tStat43_sub(:,end) = reg43_sub.tstat;   
%
R2_1_sub(end) = reg1_sub.rbar;
R2_1_sub_CPo(end) = reg1_sub_CPo.rbar;
R2_2_sub(end) = reg2_sub.rbar;
R2_2_sub_CPo(end) = reg2_sub_CPo.rbar;
R2_3_sub(end) = reg3_sub.rbar;
R2_41_sub(end) = reg41_sub.rbar;
R2_41_sub_CPo(end) = reg41_sub_CPo.rbar;
R2_42_sub(end) = reg42_sub.rbar;
R2_42_sub_CPo(end) = reg42_sub_CPo.rbar;
R2_43_sub(end) = reg43_sub.rbar;   
%
clearvars reg*

if strcmp(p.M,'ROT'), p.M=ceil((T*p.N)^(2/5)); end

% Initialize Matrices
bsBetaHat1_CG = NaN(p.B,height(betaHat1),15);
bsBetaHat1_CPo_CG = NaN(p.B,height(betaHat1),15);
bsBetaHat2_CG = NaN(p.B,height(betaHat2),15);
bsBetaHat2_CPo_CG = NaN(p.B,height(betaHat2),15);
bsBetaHat3_CG = NaN(p.B,height(betaHat3),15);
bsBetaSEHat1_CG = NaN(p.B,height(betaHat1),15);
bsBetaSEHat1_CPo_CG = NaN(p.B,height(betaHat1),15);
bsBetaSEHat2_CG = NaN(p.B,height(betaHat2),15);
bsBetaSEHat2_CPo_CG = NaN(p.B,height(betaHat2),15);
bsBetaSEHat3_CG = NaN(p.B,height(betaHat3),15);
bsBetaHat41_CG = NaN(p.B,height(betaHat41),15);
bsBetaHat41_CPo_CG = NaN(p.B,height(betaHat41),15);
bsBetaHat42_CG = NaN(p.B,height(betaHat42),15);
bsBetaHat42_CPo_CG = NaN(p.B,height(betaHat42),15);
bsBetaHat43_CG = NaN(p.B,height(betaHat43),15);
bsBetaSEHat41_CG = NaN(p.B,height(betaHat41),15);
bsBetaSEHat41_CPo_CG = NaN(p.B,height(betaHat41),15);
bsBetaSEHat42_CG = NaN(p.B,height(betaHat42),15);
bsBetaSEHat42_CPo_CG = NaN(p.B,height(betaHat42),15);
bsBetaSEHat43_CG = NaN(p.B,height(betaHat43),15);
bsR2_1_CG = NaN(p.B,15);
bsR2_1_CPo_CG = NaN(p.B,15);
bsR2_2_CG = NaN(p.B,15);
bsR2_2_CPo_CG = NaN(p.B,15);
bsR2_3_CG = NaN(p.B,15);
bsR2_41_CG = NaN(p.B,15);
bsR2_41_CPo_CG = NaN(p.B,15);
bsR2_42_CG = NaN(p.B,15);
bsR2_42_CPo_CG = NaN(p.B,15);
bsR2_43_CG = NaN(p.B,15);
%
bsBetaHat1_sub_CG = NaN(p.B,height(betaHat1_sub),15);
bsBetaHat1_sub_CPo_CG = NaN(p.B,height(betaHat1_sub),15);
bsBetaHat2_sub_CG = NaN(p.B,height(betaHat2_sub),15);
bsBetaHat2_sub_CPo_CG = NaN(p.B,height(betaHat2_sub),15);
bsBetaHat3_sub_CG = NaN(p.B,height(betaHat3_sub),15);
bsBetaSEHat1_sub_CG = NaN(p.B,height(betaHat1_sub),15);
bsBetaSEHat1_sub_CPo_CG = NaN(p.B,height(betaHat1_sub),15);
bsBetaSEHat2_sub_CG = NaN(p.B,height(betaHat2_sub),15);
bsBetaSEHat2_sub_CPo_CG = NaN(p.B,height(betaHat2_sub),15);
bsBetaSEHat3_sub_CG = NaN(p.B,height(betaHat3_sub),15);
bsBetaHat41_sub_CG = NaN(p.B,height(betaHat41_sub),15);
bsBetaHat41_sub_CPo_CG = NaN(p.B,height(betaHat41_sub),15);
bsBetaHat42_sub_CG = NaN(p.B,height(betaHat42_sub),15);
bsBetaHat42_sub_CPo_CG = NaN(p.B,height(betaHat42_sub),15);
bsBetaHat43_sub_CG = NaN(p.B,height(betaHat43_sub),15);
bsBetaSEHat41_sub_CG = NaN(p.B,height(betaHat41_sub),15);
bsBetaSEHat41_sub_CPo_CG = NaN(p.B,height(betaHat41_sub),15);
bsBetaSEHat42_sub_CG = NaN(p.B,height(betaHat42_sub),15);
bsBetaSEHat42_sub_CPo_CG = NaN(p.B,height(betaHat42_sub),15);
bsBetaSEHat43_sub_CG = NaN(p.B,height(betaHat43_sub),15);
bsR2_1_sub_CG = NaN(p.B,15);
bsR2_1_sub_CPo_CG = NaN(p.B,15);
bsR2_2_sub_CG = NaN(p.B,15);
bsR2_2_sub_CPo_CG = NaN(p.B,15);
bsR2_3_sub_CG = NaN(p.B,15);
bsR2_41_sub_CG = NaN(p.B,15);
bsR2_41_sub_CPo_CG = NaN(p.B,15);
bsR2_42_sub_CG = NaN(p.B,15);
bsR2_42_sub_CPo_CG = NaN(p.B,15);
bsR2_43_sub_CG = NaN(p.B,15);

data.fwd = NaN(size(data.yld));
data.fwd(:,1) = data.yld(:,1);
for i = 2:size(data.yld,2)
    data.fwd(:,i) = (p.yldMats(i)*data.yld(:,i) - p.yldMats(i-1)*data.yld(:,i-1))/(p.yldMats(i)-p.yldMats(i-1));
end

yldSample = data.yld(idxSampleStart:idxSampleEnd,:);
fwdSample = data.fwd(idxSampleStart:idxSampleEnd,:);

% Initial conditions for forwards (1961Q4) 
fwdSample0 = data.fwd(5,:);
% Initial conditions for yoy PCE (1961Q4) 
pred0_corePCEyoy = data.corePCEyoy(5,:);
pred0_coreCPIyoy = data.coreCPIyoy(5,:);
% Initial condition for PTR is taken from December 1961 Livingston Survey:
% Use 6m6m forward inflation at an annual rate
pred0_PTR = 100*((130.21/129.30)^2-1);

T_all = idxSampleEnd-idxSampleStart+1+40;

for b = 1:p.B

    % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % CG Bootstrap: BR Specification (PTR, r-star)
    % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%        
        % BR Specification
        bsCG1 = bootstrapCG(fwdSample,PTRSample,fwdSample0,pred0_PTR,T_all,idxSampleStart-4,idxSampleEnd-4,p);
        bsCG2 = bootstrapCG(fwdSample,[PTRSample data.corePCEyoy(idxSampleStart:idxSampleEnd)],fwdSample0,[pred0_PTR,pred0_corePCEyoy],T_all,idxSampleStart-4,idxSampleEnd-4,p);
        bsCG3 = bootstrapCG(fwdSample,data.corePCEyoy(idxSampleStart:idxSampleEnd),fwdSample0,pred0_corePCEyoy,T_all,idxSampleStart-4,idxSampleEnd-4,p);
        % CPo Specification
        bsCG1_CPo = bootstrapCG(fwdSample,data.coreCPIyoy(idxSampleStart:idxSampleEnd),fwdSample0,pred0_coreCPIyoy,T_all,idxSampleStart-4,idxSampleEnd-4,p);
        if ~o.rstarBasedOnCPI
            bsCG2_CPo = bootstrapCG(fwdSample,[data.coreCPIyoy(idxSampleStart:idxSampleEnd) data.corePCEyoy(idxSampleStart:idxSampleEnd)],fwdSample0,[pred0_coreCPIyoy,pred0_corePCEyoy],T_all,idxSampleStart-4,idxSampleEnd-4,p);
        else
            bsCG2_CPo = bootstrapCG(fwdSample,data.coreCPIyoy(idxSampleStart:idxSampleEnd),fwdSample0,pred0_coreCPIyoy,T_all,idxSampleStart-4,idxSampleEnd-4,p);
        end

        % Make r-star
        bs2Rstar  = ewma(bsCG2.bsYldCG(:,1)-bsCG2.bsPred(:,2),0.98);
        bs3Rstar  = ewma(bsCG3.bsYldCG(:,1)-bsCG3.bsPred(:,1),0.98);
        bs2Rstar_CPo = ewma(bsCG2_CPo.bsYldCG(:,1)-bsCG2_CPo.bsPred(:,end),0.98);
        % Subtract 4 because start date for datesAll is 1960Q4 whereas our bootstrap starts in "1961Q4"
        bs2RstarSample = bs2Rstar(idxSampleStart-4:idxSampleEnd-4);
        bs3RstarSample = bs3Rstar(idxSampleStart-4:idxSampleEnd-4);
        bs2RstarSample_CPo = bs2Rstar_CPo(idxSampleStart-4:idxSampleEnd-4);
        %
        bs1PTRSample = bsCG1.bsPred(idxSampleStart-4:idxSampleEnd-4,1);
        bs2PTRSample = bsCG2.bsPred(idxSampleStart-4:idxSampleEnd-4,1);
        %
        % Make CPo Factor
        bs1CPoQ  = ewma_trunc(bsCG1_CPo.bsPred(:,1),0.9615,40);
        bs2CPoQ  = ewma_trunc(bsCG2_CPo.bsPred(:,1),0.9615,40);
        bs1CPoQSample = bs1CPoQ(idxSampleStart-4:idxSampleEnd-4,1);        
        bs2CPoQSample = bs2CPoQ(idxSampleStart-4:idxSampleEnd-4,1);                
        %
        bs1XrnSample = bsCG1.bsXrn(idxSampleStart-4:idxSampleEnd-4,4:end);
        bs1Xrn4Sample = bsCG1.bsXrn4(idxSampleStart-4:idxSampleEnd-4,4:end);
        bs2XrnSample = bsCG2.bsXrn(idxSampleStart-4:idxSampleEnd-4,4:end);
        bs2Xrn4Sample = bsCG2.bsXrn4(idxSampleStart-4:idxSampleEnd-4,4:end);
        bs3XrnSample = bsCG3.bsXrn(idxSampleStart-4:idxSampleEnd-4,4:end);
        bs3Xrn4Sample = bsCG3.bsXrn4(idxSampleStart-4:idxSampleEnd-4,4:end);
        bs1XrnSample_CPo = bsCG1_CPo.bsXrn(idxSampleStart-4:idxSampleEnd-4,4:end);
        bs1Xrn4Sample_CPo = bsCG1_CPo.bsXrn4(idxSampleStart-4:idxSampleEnd-4,4:end);
        bs2XrnSample_CPo = bsCG2_CPo.bsXrn(idxSampleStart-4:idxSampleEnd-4,4:end);
        bs2Xrn4Sample_CPo = bsCG2_CPo.bsXrn4(idxSampleStart-4:idxSampleEnd-4,4:end);
        %
        bs1PCsscSample = bsCG1.bsPCssc(idxSampleStart-4:idxSampleEnd-4,:);
        bs2PCsscSample = bsCG2.bsPCssc(idxSampleStart-4:idxSampleEnd-4,:);
        bs3PCsscSample = bsCG3.bsPCssc(idxSampleStart-4:idxSampleEnd-4,:);
        bs1PCsscSample_CPo = bsCG1_CPo.bsPCssc(idxSampleStart-4:idxSampleEnd-4,:);
        bs2PCsscSample_CPo = bsCG2_CPo.bsPCssc(idxSampleStart-4:idxSampleEnd-4,:);

        for i = 1:14

            bsReg1 = nwest(bs1XrnSample(2:end,i),[ones(T-1,1), bs1PCsscSample(1:end-1,:) bs1PTRSample(1:end-1)],0);
            bsReg41 = nwest(bs1Xrn4Sample(5:end,i),[ones(T-4,1), bs1PCsscSample(1:end-4,:), bs1PTRSample(1:end-4)],nwLagh4); 
            bsReg2 = nwest(bs2XrnSample(2:end,i),[ones(T-1,1), bs2PCsscSample(1:end-1,:), bs2PTRSample(1:end-1) bs2RstarSample(1:end-1)],0);
            bsReg42 = nwest(bs2Xrn4Sample(5:end,i),[ones(T-4,1), bs2PCsscSample(1:end-4,:) bs2PTRSample(1:end-4) bs2RstarSample(1:end-4)],nwLagh4); 
            bsReg3 = nwest(bs3XrnSample(2:end,i),[ones(T-1,1), bs3PCsscSample(1:end-1,:), bs3RstarSample(1:end-1)],0);            
            bsReg43 = nwest(bs3Xrn4Sample(5:end,i),[ones(T-4,1), bs3PCsscSample(1:end-4,:) bs3RstarSample(1:end-4)],nwLagh4);
            %
            bsReg1_CPo = nwest(bs1XrnSample_CPo(2:end,i),[ones(T-1,1), bs1PCsscSample_CPo(1:end-1,:) bs1CPoQSample(1:end-1)],0);
            bsReg41_CPo = nwest(bs1Xrn4Sample_CPo(5:end,i),[ones(T-4,1), bs1PCsscSample_CPo(1:end-4,:), bs1CPoQSample(1:end-4)],nwLagh4); 
            bsReg2_CPo = nwest(bs2XrnSample_CPo(2:end,i),[ones(T-1,1), bs2PCsscSample_CPo(1:end-1,:) bs2CPoQSample(1:end-1) bs2RstarSample_CPo(1:end-1)],0);
            bsReg42_CPo = nwest(bs2Xrn4Sample_CPo(5:end,i),[ones(T-4,1), bs2PCsscSample_CPo(1:end-4,:), bs2CPoQSample(1:end-4) bs2RstarSample_CPo(1:end-4)],nwLagh4);         
    
            bsBetaHat1_CG(b,:,i) = bsReg1.beta;
            bsBetaHat2_CG(b,:,i) = bsReg2.beta;
            bsBetaHat3_CG(b,:,i) = bsReg3.beta;
            bsBetaSEHat1_CG(b,:,i) = bsReg1.beta./bsReg1.tstat;
            bsBetaSEHat2_CG(b,:,i) = bsReg2.beta./bsReg2.tstat;
            bsBetaSEHat3_CG(b,:,i) = bsReg3.beta./bsReg3.tstat;
            bsBetaHat1_CPo_CG(b,:,i) = bsReg1_CPo.beta;
            bsBetaHat2_CPo_CG(b,:,i) = bsReg2_CPo.beta;
            bsBetaSEHat1_CPo_CG(b,:,i) = bsReg1_CPo.beta./bsReg1_CPo.tstat;
            bsBetaSEHat2_CPo_CG(b,:,i) = bsReg2_CPo.beta./bsReg2_CPo.tstat;        
            %
            bsBetaHat41_CG(b,:,i) = bsReg41.beta;
            bsBetaHat42_CG(b,:,i) = bsReg42.beta;
            bsBetaHat43_CG(b,:,i) = bsReg43.beta;
            bsBetaSEHat41_CG(b,:,i) = bsReg41.beta./bsReg41.tstat;
            bsBetaSEHat42_CG(b,:,i) = bsReg42.beta./bsReg42.tstat;
            bsBetaSEHat43_CG(b,:,i) = bsReg43.beta./bsReg43.tstat;
            bsBetaHat41_CPo_CG(b,:,i) = bsReg41_CPo.beta;
            bsBetaHat42_CPo_CG(b,:,i) = bsReg42_CPo.beta;
            bsBetaSEHat41_CPo_CG(b,:,i) = bsReg41_CPo.beta./bsReg41_CPo.tstat;
            bsBetaSEHat42_CPo_CG(b,:,i) = bsReg42_CPo.beta./bsReg42_CPo.tstat;      
            %
            bsR2_1_CG(b,i) = bsReg1.rbar;
            bsR2_1_CPo_CG(b,i) = bsReg1_CPo.rbar;
            bsR2_2_CG(b,i) = bsReg2.rbar;
            bsR2_2_CPo_CG(b,i) = bsReg2_CPo.rbar;
            bsR2_3_CG(b,i) = bsReg3.rbar;
            bsR2_41_CG(b,i) = bsReg41.rbar;
            bsR2_41_CPo_CG(b,i) = bsReg41_CPo.rbar;
            bsR2_42_CG(b,i) = bsReg42.rbar;
            bsR2_42_CPo_CG(b,i) = bsReg42_CPo.rbar;
            bsR2_43_CG(b,i) = bsReg43.rbar;      
            %
            clearvars reg*
            %
            bsReg1_sub = nwest(bs1XrnSample(2+offset:end,i),[ones(T-1-offset,1), bs1PCsscSample(1+offset:end-1,:) bs1PTRSample(1+offset:end-1)],0);
            bsReg41_sub = nwest(bs1Xrn4Sample(5+offset:end,i),[ones(T-4-offset,1), bs1PCsscSample(1+offset:end-4,:), bs1PTRSample(1+offset:end-4)],nwLagh4); 
            bsReg2_sub = nwest(bs2XrnSample(2+offset:end,i),[ones(T-1-offset,1), bs2PCsscSample(1+offset:end-1,:), bs2PTRSample(1+offset:end-1) bs2RstarSample(1+offset:end-1)],0);
            bsReg42_sub = nwest(bs2Xrn4Sample(5+offset:end,i),[ones(T-4-offset,1), bs2PCsscSample(1+offset:end-4,:) bs2PTRSample(1+offset:end-4) bs2RstarSample(1+offset:end-4)],nwLagh4); 
            bsReg3_sub = nwest(bs3XrnSample(2+offset:end,i),[ones(T-1-offset,1), bs3PCsscSample(1+offset:end-1,:), bs3RstarSample(1+offset:end-1)],0);
            bsReg43_sub = nwest(bs3Xrn4Sample(5+offset:end,i),[ones(T-4-offset,1), bs3PCsscSample(1+offset:end-4,:) bs3RstarSample(1+offset:end-4)],nwLagh4);
            %
            bsReg1_sub_CPo = nwest(bs1XrnSample_CPo(2+offset:end,i),[ones(T-1-offset,1), bs1PCsscSample_CPo(1+offset:end-1,:) bs1CPoQSample(1+offset:end-1)],0);
            bsReg41_sub_CPo = nwest(bs1Xrn4Sample_CPo(5+offset:end,i),[ones(T-4-offset,1), bs1PCsscSample_CPo(1+offset:end-4,:), bs1CPoQSample(1+offset:end-4)],nwLagh4); 
            bsReg2_sub_CPo = nwest(bs2XrnSample_CPo(2+offset:end,i),[ones(T-1-offset,1), bs2PCsscSample_CPo(1+offset:end-1,:) bs2CPoQSample(1+offset:end-1) bs2RstarSample_CPo(1+offset:end-1)],0);
            bsReg42_sub_CPo = nwest(bs2Xrn4Sample_CPo(5+offset:end,i),[ones(T-4-offset,1), bs2PCsscSample_CPo(1+offset:end-4,:), bs2CPoQSample(1+offset:end-4) bs2RstarSample_CPo(1+offset:end-4)],nwLagh4);         
            %
            bsBetaHat1_sub_CG(b,:,i) = bsReg1_sub.beta;
            bsBetaHat2_sub_CG(b,:,i) = bsReg2_sub.beta;
            bsBetaHat3_sub_CG(b,:,i) = bsReg3_sub.beta;
            bsBetaSEHat1_sub_CG(b,:,i) = bsReg1_sub.beta./bsReg1_sub.tstat;
            bsBetaSEHat2_sub_CG(b,:,i) = bsReg2_sub.beta./bsReg2_sub.tstat;
            bsBetaSEHat3_sub_CG(b,:,i) = bsReg3_sub.beta./bsReg3_sub.tstat;
            bsBetaHat1_sub_CPo_CG(b,:,i) = bsReg1_sub_CPo.beta;
            bsBetaHat2_sub_CPo_CG(b,:,i) = bsReg2_sub_CPo.beta;
            bsBetaSEHat1_sub_CPo_CG(b,:,i) = bsReg1_sub_CPo.beta./bsReg1_sub_CPo.tstat;
            bsBetaSEHat2_sub_CPo_CG(b,:,i) = bsReg2_sub_CPo.beta./bsReg2_sub_CPo.tstat;        
            %
            bsBetaHat41_sub_CG(b,:,i) = bsReg41_sub.beta;
            bsBetaHat42_sub_CG(b,:,i) = bsReg42_sub.beta;
            bsBetaHat43_sub_CG(b,:,i) = bsReg43_sub.beta;
            bsBetaSEHat41_sub_CG(b,:,i) = bsReg41_sub.beta./bsReg41_sub.tstat;
            bsBetaSEHat42_sub_CG(b,:,i) = bsReg42_sub.beta./bsReg42_sub.tstat;
            bsBetaSEHat43_sub_CG(b,:,i) = bsReg43_sub.beta./bsReg43_sub.tstat;
            bsBetaHat41_sub_CPo_CG(b,:,i) = bsReg41_sub_CPo.beta;
            bsBetaHat42_sub_CPo_CG(b,:,i) = bsReg42_sub_CPo.beta;
            bsBetaSEHat41_sub_CPo_CG(b,:,i) = bsReg41_sub_CPo.beta./bsReg41_sub_CPo.tstat;
            bsBetaSEHat42_sub_CPo_CG(b,:,i) = bsReg42_sub_CPo.beta./bsReg42_sub_CPo.tstat;      
            %
            bsR2_1_sub_CG(b,i) = bsReg1_sub.rbar;
            bsR2_1_sub_CPo_CG(b,i) = bsReg1_sub_CPo.rbar;
            bsR2_2_sub_CG(b,i) = bsReg2_sub.rbar;
            bsR2_2_sub_CPo_CG(b,i) = bsReg2_sub_CPo.rbar;
            bsR2_3_sub_CG(b,i) = bsReg3_sub.rbar;
            bsR2_41_sub_CG(b,i) = bsReg41_sub.rbar;
            bsR2_41_sub_CPo_CG(b,i) = bsReg41_sub_CPo.rbar;
            bsR2_42_sub_CG(b,i) = bsReg42_sub.rbar;
            bsR2_42_sub_CPo_CG(b,i) = bsReg42_sub_CPo.rbar;
            bsR2_43_sub_CG(b,i) = bsReg43_sub.rbar;   
            %
            clearvars reg*
        end

        bsReg1 = nwest(mean(bs1XrnSample(2:end,:),2),[ones(T-1,1), bs1PCsscSample(1:end-1,:) bs1PTRSample(1:end-1)],0);
        bsReg41 = nwest(mean(bs1Xrn4Sample(5:end,:),2),[ones(T-4,1), bs1PCsscSample(1:end-4,:), bs1PTRSample(1:end-4)],nwLagh4); 
        bsReg2 = nwest(mean(bs2XrnSample(2:end,:),2),[ones(T-1,1), bs2PCsscSample(1:end-1,:), bs2PTRSample(1:end-1) bs2RstarSample(1:end-1)],0);
        bsReg42 = nwest(mean(bs2Xrn4Sample(5:end,:),2),[ones(T-4,1), bs2PCsscSample(1:end-4,:) bs2PTRSample(1:end-4) bs2RstarSample(1:end-4)],nwLagh4); 
        bsReg3 = nwest(mean(bs3XrnSample(2:end,:),2),[ones(T-1,1), bs3PCsscSample(1:end-1,:), bs3RstarSample(1:end-1)],0);
        bsReg43 = nwest(mean(bs3Xrn4Sample(5:end,:),2),[ones(T-4,1), bs3PCsscSample(1:end-4,:) bs3RstarSample(1:end-4)],nwLagh4);
        bsReg1_CPo = nwest(mean(bs1XrnSample_CPo(2:end,:),2),[ones(T-1,1), bs1PCsscSample_CPo(1:end-1,:) bs1CPoQSample(1:end-1)],0);
        bsReg41_CPo = nwest(mean(bs1Xrn4Sample_CPo(5:end,:),2),[ones(T-4,1), bs1PCsscSample_CPo(1:end-4,:), bs1CPoQSample(1:end-4)],nwLagh4); 
        bsReg2_CPo = nwest(mean(bs2XrnSample_CPo(2:end,:),2),[ones(T-1,1), bs2PCsscSample_CPo(1:end-1,:) bs2CPoQSample(1:end-1) bs2RstarSample_CPo(1:end-1)],0);
        bsReg42_CPo = nwest(mean(bs2Xrn4Sample_CPo(5:end,:),2),[ones(T-4,1), bs2PCsscSample_CPo(1:end-4,:), bs2CPoQSample(1:end-4) bs2RstarSample_CPo(1:end-4)],nwLagh4);         

        bsBetaHat1_CG(b,:,end) = bsReg1.beta;
        bsBetaHat2_CG(b,:,end) = bsReg2.beta;
        bsBetaHat3_CG(b,:,end) = bsReg3.beta;
        bsBetaSEHat1_CG(b,:,end) = bsReg1.beta./bsReg1.tstat;
        bsBetaSEHat2_CG(b,:,end) = bsReg2.beta./bsReg2.tstat;
        bsBetaSEHat3_CG(b,:,end) = bsReg3.beta./bsReg3.tstat;
        bsBetaHat1_CPo_CG(b,:,end) = bsReg1_CPo.beta;
        bsBetaHat2_CPo_CG(b,:,end) = bsReg2_CPo.beta;
        bsBetaSEHat1_CPo_CG(b,:,end) = bsReg1_CPo.beta./bsReg1_CPo.tstat;
        bsBetaSEHat2_CPo_CG(b,:,end) = bsReg2_CPo.beta./bsReg2_CPo.tstat;        
        %
        bsBetaHat41_CG(b,:,end) = bsReg41.beta;
        bsBetaHat42_CG(b,:,end) = bsReg42.beta;
        bsBetaHat43_CG(b,:,end) = bsReg43.beta;
        bsBetaSEHat41_CG(b,:,end) = bsReg41.beta./bsReg41.tstat;
        bsBetaSEHat42_CG(b,:,end) = bsReg42.beta./bsReg42.tstat;
        bsBetaSEHat43_CG(b,:,end) = bsReg43.beta./bsReg43.tstat;
        bsBetaHat41_CPo_CG(b,:,end) = bsReg41_CPo.beta;
        bsBetaHat42_CPo_CG(b,:,end) = bsReg42_CPo.beta;
        bsBetaSEHat41_CPo_CG(b,:,end) = bsReg41_CPo.beta./bsReg41_CPo.tstat;
        bsBetaSEHat42_CPo_CG(b,:,end) = bsReg42_CPo.beta./bsReg42_CPo.tstat;        
        %
        bsR2_1_CG(b,end) = bsReg1.rbar;
        bsR2_2_CG(b,end) = bsReg2.rbar;
        bsR2_3_CG(b,end) = bsReg3.rbar;
        bsR2_41_CG(b,end) = bsReg41.rbar;
        bsR2_42_CG(b,end) = bsReg42.rbar;
        bsR2_43_CG(b,end) = bsReg43.rbar;
        bsR2_1_CPo_CG(b,end) = bsReg1_CPo.rbar;
        bsR2_2_CPo_CG(b,end) = bsReg2_CPo.rbar;
        bsR2_41_CPo_CG(b,end) = bsReg41_CPo.rbar;
        bsR2_42_CPo_CG(b,end) = bsReg42_CPo.rbar;  
        %
        clearvars reg*
        %
        bsReg1_sub = nwest(mean(bs1XrnSample(2+offset:end,:),2),[ones(T-1-offset,1), bs1PCsscSample(1+offset:end-1,:) bs1PTRSample(1+offset:end-1)],0);
        bsReg41_sub = nwest(mean(bs1Xrn4Sample(5+offset:end,:),2),[ones(T-4-offset,1), bs1PCsscSample(1+offset:end-4,:), bs1PTRSample(1+offset:end-4)],nwLagh4); 
        bsReg2_sub = nwest(mean(bs2XrnSample(2+offset:end,:),2),[ones(T-1-offset,1), bs2PCsscSample(1+offset:end-1,:), bs2PTRSample(1+offset:end-1) bs2RstarSample(1+offset:end-1)],0);
        bsReg42_sub = nwest(mean(bs2Xrn4Sample(5+offset:end,:),2),[ones(T-4-offset,1), bs2PCsscSample(1+offset:end-4,:) bs2PTRSample(1+offset:end-4) bs2RstarSample(1+offset:end-4)],nwLagh4); 
        bsReg3_sub = nwest(mean(bs3XrnSample(2+offset:end,:),2),[ones(T-1-offset,1), bs3PCsscSample(1+offset:end-1,:), bs3RstarSample(1+offset:end-1)],0);
        bsReg43_sub = nwest(mean(bs3Xrn4Sample(5+offset:end,:),2),[ones(T-4-offset,1), bs3PCsscSample(1+offset:end-4,:) bs3RstarSample(1+offset:end-4)],nwLagh4);
        bsReg1_sub_CPo = nwest(mean(bs1XrnSample_CPo(2+offset:end,:),2),[ones(T-1-offset,1), bs1PCsscSample_CPo(1+offset:end-1,:) bs1CPoQSample(1+offset:end-1)],0);
        bsReg41_sub_CPo = nwest(mean(bs1Xrn4Sample_CPo(5+offset:end,:),2),[ones(T-4-offset,1), bs1PCsscSample_CPo(1+offset:end-4,:), bs1CPoQSample(1+offset:end-4)],nwLagh4); 
        bsReg2_sub_CPo = nwest(mean(bs2XrnSample_CPo(2+offset:end,:),2),[ones(T-1-offset,1), bs2PCsscSample_CPo(1+offset:end-1,:) bs2CPoQSample(1+offset:end-1) bs2RstarSample_CPo(1+offset:end-1)],0);
        bsReg42_sub_CPo = nwest(mean(bs2Xrn4Sample_CPo(5+offset:end,:),2),[ones(T-4-offset,1), bs2PCsscSample_CPo(1+offset:end-4,:), bs2CPoQSample(1+offset:end-4) bs2RstarSample_CPo(1+offset:end-4)],nwLagh4);         
        %
        bsBetaHat1_sub_CG(b,:,end) = bsReg1_sub.beta;
        bsBetaHat2_sub_CG(b,:,end) = bsReg2_sub.beta;
        bsBetaHat3_sub_CG(b,:,end) = bsReg3_sub.beta;
        bsBetaSEHat1_sub_CG(b,:,end) = bsReg1_sub.beta./bsReg1_sub.tstat;
        bsBetaSEHat2_sub_CG(b,:,end) = bsReg2_sub.beta./bsReg2_sub.tstat;
        bsBetaSEHat3_sub_CG(b,:,end) = bsReg3_sub.beta./bsReg3_sub.tstat;
        bsBetaHat1_sub_CPo_CG(b,:,end) = bsReg1_sub_CPo.beta;
        bsBetaHat2_sub_CPo_CG(b,:,end) = bsReg2_sub_CPo.beta;
        bsBetaSEHat1_sub_CPo_CG(b,:,end) = bsReg1_sub_CPo.beta./bsReg1_sub_CPo.tstat;
        bsBetaSEHat2_sub_CPo_CG(b,:,end) = bsReg2_sub_CPo.beta./bsReg2_sub_CPo.tstat;        
        %
        bsBetaHat41_sub_CG(b,:,end) = bsReg41_sub.beta;
        bsBetaHat42_sub_CG(b,:,end) = bsReg42_sub.beta;
        bsBetaHat43_sub_CG(b,:,end) = bsReg43_sub.beta;
        bsBetaSEHat41_sub_CG(b,:,end) = bsReg41_sub.beta./bsReg41_sub.tstat;
        bsBetaSEHat42_sub_CG(b,:,end) = bsReg42_sub.beta./bsReg42_sub.tstat;
        bsBetaSEHat43_sub_CG(b,:,end) = bsReg43_sub.beta./bsReg43_sub.tstat;
        bsBetaHat41_sub_CPo_CG(b,:,end) = bsReg41_sub_CPo.beta;
        bsBetaHat42_sub_CPo_CG(b,:,end) = bsReg42_sub_CPo.beta;
        bsBetaSEHat41_sub_CPo_CG(b,:,end) = bsReg41_sub_CPo.beta./bsReg41_sub_CPo.tstat;
        bsBetaSEHat42_sub_CPo_CG(b,:,end) = bsReg42_sub_CPo.beta./bsReg42_sub_CPo.tstat;        
        %
        bsR2_1_sub_CG(b,end) = bsReg1_sub.rbar;
        bsR2_2_sub_CG(b,end) = bsReg2_sub.rbar;
        bsR2_3_sub_CG(b,end) = bsReg3_sub.rbar;
        bsR2_41_sub_CG(b,end) = bsReg41_sub.rbar;
        bsR2_42_sub_CG(b,end) = bsReg42_sub.rbar;
        bsR2_43_sub_CG(b,end) = bsReg43_sub.rbar;
        bsR2_1_sub_CPo_CG(b,end) = bsReg1_sub_CPo.rbar;
        bsR2_2_sub_CPo_CG(b,end) = bsReg2_sub_CPo.rbar;
        bsR2_41_sub_CPo_CG(b,end) = bsReg41_sub_CPo.rbar;
        bsR2_42_sub_CPo_CG(b,end) = bsReg42_sub_CPo.rbar;
        %
        clearvars reg*     

        if mod(b,50)==0, disp([num2str(b) ' of ' num2str(p.B), ' completed']); end

end

pval1 = NaN(size(betaHat1));
pval2 = NaN(size(betaHat2));
pval3 = NaN(size(betaHat3));
pval1_CPo = NaN(size(betaHat1));
pval2_CPo = NaN(size(betaHat2));
pval41 = NaN(size(betaHat1));
pval42 = NaN(size(betaHat2));
pval43 = NaN(size(betaHat3));
pval41_CPo = NaN(size(betaHat1));
pval42_CPo = NaN(size(betaHat2));
%
pval1_sub = NaN(size(betaHat1));
pval2_sub = NaN(size(betaHat2));
pval3_sub = NaN(size(betaHat3));
pval1_sub_CPo = NaN(size(betaHat1));
pval2_sub_CPo = NaN(size(betaHat2));
pval41_sub = NaN(size(betaHat1));
pval42_sub = NaN(size(betaHat2));
pval43_sub = NaN(size(betaHat3));
pval41_sub_CPo = NaN(size(betaHat1));
pval42_sub_CPo = NaN(size(betaHat2));

for i = 1:15
    pval1(:,i) = makeBsPval(bsBetaHat1_CG(:,:,i),betaHat1(:,i),bsBetaSEHat1_CG(:,:,i),tStat1(:,i),1);
    pval2(:,i) = makeBsPval(bsBetaHat2_CG(:,:,i),betaHat2(:,i),bsBetaSEHat2_CG(:,:,i),tStat2(:,i),1);
    pval3(:,i) = makeBsPval(bsBetaHat3_CG(:,:,i),betaHat3(:,i),bsBetaSEHat3_CG(:,:,i),tStat3(:,i),1);
    pval1_CPo(:,i) = makeBsPval(bsBetaHat1_CPo_CG(:,:,i),betaHat1_CPo(:,i),bsBetaSEHat1_CPo_CG(:,:,i),tStat1_CPo(:,i),1);
    pval2_CPo(:,i) = makeBsPval(bsBetaHat2_CPo_CG(:,:,i),betaHat2_CPo(:,i),bsBetaSEHat2_CPo_CG(:,:,i),tStat2_CPo(:,i),1);
    %
    pval41(:,i) = makeBsPval(bsBetaHat41_CG(:,:,i),betaHat41(:,i),bsBetaSEHat41_CG(:,:,i),tStat41(:,i),1);
    pval42(:,i) = makeBsPval(bsBetaHat42_CG(:,:,i),betaHat42(:,i),bsBetaSEHat42_CG(:,:,i),tStat42(:,i),1);
    pval43(:,i) = makeBsPval(bsBetaHat43_CG(:,:,i),betaHat43(:,i),bsBetaSEHat43_CG(:,:,i),tStat43(:,i),1);
    pval41_CPo(:,i) = makeBsPval(bsBetaHat41_CPo_CG(:,:,i),betaHat41_CPo(:,i),bsBetaSEHat41_CPo_CG(:,:,i),tStat41_CPo(:,i),1);
    pval42_CPo(:,i) = makeBsPval(bsBetaHat42_CPo_CG(:,:,i),betaHat42_CPo(:,i),bsBetaSEHat42_CPo_CG(:,:,i),tStat42_CPo(:,i),1);
    %
    pval1_sub(:,i) = makeBsPval(bsBetaHat1_sub_CG(:,:,i),betaHat1_sub(:,i),bsBetaSEHat1_sub_CG(:,:,i),tStat1_sub(:,i),1);
    pval2_sub(:,i) = makeBsPval(bsBetaHat2_sub_CG(:,:,i),betaHat2_sub(:,i),bsBetaSEHat2_sub_CG(:,:,i),tStat2_sub(:,i),1);
    pval3_sub(:,i) = makeBsPval(bsBetaHat3_sub_CG(:,:,i),betaHat3_sub(:,i),bsBetaSEHat3_sub_CG(:,:,i),tStat3_sub(:,i),1);
    pval1_sub_CPo(:,i) = makeBsPval(bsBetaHat1_sub_CPo_CG(:,:,i),betaHat1_sub_CPo(:,i),bsBetaSEHat1_sub_CPo_CG(:,:,i),tStat1_sub_CPo(:,i),1);
    pval2_sub_CPo(:,i) = makeBsPval(bsBetaHat2_sub_CPo_CG(:,:,i),betaHat2_sub_CPo(:,i),bsBetaSEHat2_sub_CPo_CG(:,:,i),tStat2_sub_CPo(:,i),1);
    %
    pval41_sub(:,i) = makeBsPval(bsBetaHat41_sub_CG(:,:,i),betaHat41_sub(:,i),bsBetaSEHat41_sub_CG(:,:,i),tStat41_sub(:,i),1);
    pval42_sub(:,i) = makeBsPval(bsBetaHat42_sub_CG(:,:,i),betaHat42_sub(:,i),bsBetaSEHat42_sub_CG(:,:,i),tStat42_sub(:,i),1);
    pval43_sub(:,i) = makeBsPval(bsBetaHat43_sub_CG(:,:,i),betaHat43_sub(:,i),bsBetaSEHat43_sub_CG(:,:,i),tStat43_sub(:,i),1);
    pval41_sub_CPo(:,i) = makeBsPval(bsBetaHat41_sub_CPo_CG(:,:,i),betaHat41_sub_CPo(:,i),bsBetaSEHat41_sub_CPo_CG(:,:,i),tStat41_sub_CPo(:,i),1);
    pval42_sub_CPo(:,i) = makeBsPval(bsBetaHat42_sub_CPo_CG(:,:,i),betaHat42_sub_CPo(:,i),bsBetaSEHat42_sub_CPo_CG(:,:,i),tStat42_sub_CPo(:,i),1);    
    
end

close all;

if idxSaveData
    tmpFileName = ['BR_IndMats_B',num2str(p.B) ,'_BRdata', num2str(o.exactDataBR), '_rstarCPI', num2str(o.rstarBasedOnCPI), '_', num2str(o.startDateSample),'_', num2str(o.endDateSample)];
    save(['../output/', tmpFileName, '.mat']);
end

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FUNCTIONS
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function out = ewma(y,a)
        T = numel(y);
        out = NaN(T,1);
        for t = 1:T
            if t==1
                out(t) = y(t);
            else
                out(t) = out(t-1) + (1-a)*(y(t) - out(t-1));
            end
        end
end

function x=ewma_trunc(y,alpha,fobs)
    n = length(y);
    weights = alpha.^(fobs-1:-1:0)';
    weights = weights./sum(weights); % normalize to sum up to one
    gap=zeros(n,1);
    for j=fobs:n
        gap(j) = sum(weights.*y(j-fobs+1:j));
    end
    x = gap;
end

